layer_boundaries_ for uniform stack lattice with 2 elements of pitch = 1.5:

(gdb) p layer_boundaries_
$1 = std::vector of length 3, capacity 4 = {0, 1.5, 3}

layer_boundaries_ for nonunifrom stack lattice with 2 elements of 'pitch' [1.0, 2.0]:

(gdb) p layer_boundaries_
$1 = std::vector of length 3, capacity 4 = {0, 1, 2}

The uniform case is correct, but the nonuniform case should be {0, 1, 3}. The same happens in the Python API so I'll need to adjust the implementation for creating the _layer_boundaries array there as well.

Note that this is because my current implementation does not catch the case where a non-iterable pitch is set for a non-uniform stack lattice. I'm working on a solution to this.

UPDATE: Commit 2bff4cc fixes the implementation and now catches errors when setting pitch for a StackLattice object.

I did some more poking around today and realized that the get_local_position function doesn't handle cases where the lattice index is invalid, which I speculate is a component of telling OpenMC that the particle is no longer in the lattice. Regardless, I added control flow to handle these cases in 6789365 and the particle not found error has disappeared.

This is a pretty big PR so I think 2 or 3 reviewers would be ideal @paulromano @pshriwise @anyone else who wants to review

Thanks @yardasol for this PR! I'm definitely planning on reviewing this and would appreciate getting @pshriwise's input too. I unfortunately have to plead with you to be patient. We're in the midst of preparing for our 4-day OpenMC course next week, hence the pileup of unreviewed PRs. Once things are all set for the course, we should be able to start handling a lot of the PRs, including this one.

Hi @paulromano, I apologize for the frequent pinging. I don't expect anyone to review PRs immediately when I tag them, so I'm a bit embarrassed that I came off in that way! I'll try to be more suave about that. I pinged you and Patrick because this PR has been sitting in draft mode for weeks now and I just wanted to make sure you knew it existed :). Good luck with the course next week!

@yardasol I apologize that I haven't gotten you a full review on this yet. This is quite an impressive PR! But, it's also a big one that has some design implications so I need a bit of time to work through it and also try out the feature as well so I can better understand it. I just wanted you to know I haven't forgotten it 😄 If you happen to be free tomorrow at 9am CT and are up for it, it would be great if you could join our OpenMC monthly meet-up and walk through an example with the stack lattice feature.

Hi @paulromano. Yes this PR is quite large. I appreciate you taking the time to go through it. What fortuitous timing for the monthly meet-up; I'd be happy to swing by and go through a notebook example with everyone.

Update: So, I decided to go with bit more sophisticated of a model for the example I want to show (a toy BWR fuel rod). Turns out that this contribution isn't ready at all in the way I've implemented it. The toy BWR fuel rod model is giving me bunch of particle not found errors again. Here's my hypothesis: StackLattice tiles are not explicitly bounded in two out of three dimensions, so if a model uses StackLattice but doesn't assign a cell to every possible point in each tile, there will be errors. The reason my integration tests work is because I happened to stumble upon this exact corner case where the current implementation works.

I'm thinking about how to allow users to specify a bound on their lattice in the unbounded dimensions. Perhaps an attribute tile_boundary that contains a Surface object (in the nonuniform case it would be an array of Surface objects) defining the boundary of each tile in the uniform case? This way we could easily compute the distance to the surface (on both the C++ backend and the Python frontend) and then if that surface is crossed, we are no longer in the lattice.

Even though I realize this is now not working as intended, I think I'll still show up to the meeting tomorrow anyways so I can get ideas and feedback from folks.

Thanks for proposing this feature @yardasol! I've had a chance to look into the implementation a little bit more and do feel that the right path forward is to try to extend our current RectLattice class to support a one-dimensional lattice (as we discussed on the last monthly meeting call). The only case that won't cover is the non-uniform case, but I think that is much less common and raises some thorny design issues. Let me know what you think and if there's any support I can provide in moving toward the 1D RectLattice capability.

Hi @paulromano, I generally agree with your conclusion that extending RectLattice to have true 1D lattices would be the best path forwards.

We'll still need to figure out how to implement the 1D tiling structure without needing to define every region in space off to infinity within the pitch-defined bounds of the tile. I'll try to think about some good ways to do this.

@yardasol have you thought more about a way to implement this? I recently built some serpent inputs using the stack lattice feature with variably spaced planes, and wow, it is so much more convenient than what has to be done in openmc at the moment.

I guess it would be something like: RectLattice would look like the RectilinearMesh instead now, and save some 1D arrays of the X, Y, and Z values the planes are defined at. I think it's a computationally negligible part of the code, so the extra work to look up into those (usually fitting in the cache) vectors would be pretty quick.

Then only in the python side, we could implement ZStack, YStack, etc as shorthand ways to define a RectLattice with the expected syntax.